      SUBROUTINE ROMUL(litter_in, Z, T, WLL,
     & cpool_soil, npool_soil, cflux, dt)

      use INMSOIL_NUMPARAMS, only : 
     & ML, MS
       USE CBALANCE_CONST
      !INCLUDE 'soil.inc'
!       PARAMETER (nv2=7) !number of vegetation types
!       input
        real (KIND=8) litter_in(nv2) !total litter input kgCm-2s-1
        real (KIND=4) Z(ML) !soil layer depth, cm
        real (KIND=4) T(ML) !soil temperature degC
        real (KIND=4) WLL(ML) !soil water content (g/g)
        
          real (KIND=8) cpool_soil(nsp), npool_soil(nsp) !also output C and N 
                                                       !in soil pools kgCm-2
                                                       !output
        real (KIND=8) cflux ! CO2 flux from the soil kgCm-2s-1


!       Local variables
        integer p
        real bbact, blumb
        PARAMETER (bbact=24.0/2.0, blumb=12.8/2.0)
        real row_lit, row_org
        PARAMETER (row_lit=0.12, row_org=0.12) 
        real leaf_litr, root_litr
        PARAMETER (leaf_litr=0.5, root_litr=0.5)
        real(8)  dt
        real deltat, mt
        PARAMETER (mt=0.8)
        real (KIND=8) csoilL_leaf,csoilL_root,  csoilHs, 
     & csoilHu, csoilF_leaf, csoilF_root, csoilF, csoilL, 
     & nsoilL_leaf, nsoilL_root,  nsoilHs, nsoilHu, 
     & nsoilF_leaf,nsoilF_root, nsoilF, nsoilL
        real acc_leaf,acc_root
      real acc_n_leaf, acc_n_root
        real n_cont_leaf0, n_cont_root0, n_cont_F0
        real ash_cont_leaf0, ash_cont_root0, ash_cont_F0
      real n_cont_leaf, n_cont_root, n_cont_F, n_cont_F_leaf,
     & n_cont_F_root
        real k1_som_leaf, k3_som_leaf,k1_som_root, k3_som_root,
     & k2_som_leaf,k2_som_root,k4_som,k5_som,
     & k6_som, k7_som, k4_som_a,k4_som_b
      real Ml_leaf,Ml_root, Mf, Mf_leaf, Mf_root, Mhs, Mhu 

      real temp_lit, temp_org, temp_min
        real masswc_lit, masswc_org, masswc_min

      real litter_inp_leaf,litter_inp_root
      real cton_leaf(nv2),cton_root(nv2),cton_sap(nv2)
      real depth_b 
      PARAMETER (depth_b=70.0) 
      real depth_int
       deltat=dt/(24.*3600.) !time step in days


       temp_lit=T(5) !surface
       temp_org=T(8)  !8 cm
       masswc_lit=wll(5)*100*1.2/row_lit 
       masswc_org=wll(8)*100*1.2/row_org 
       k=ms+1
       masswc_min=0.0
       temp_min=0.0
       k6_som=0.0
       depth_int=0.0
       do while (depth_b.gt.z(k))
       masswc_min=masswc_min+0.5*wll(k)*(z(k+1)-z(k-1))
       temp_min=temp_min+0.5*T(k)*(z(k+1)-z(k-1))
!         k6_som=k6_som+0.5*KHUM(k)*(z(k+1)-z(k-1))
       depth_int=depth_int+0.5*(z(k+1)-z(k-1))
       k=k+1
       end do
       masswc_min=100*masswc_min/depth_int  
       temp_min=temp_min/depth_int  
!         k6_som=k6_som/depth_int
       k6_som=0.00006
!**************************initialize
       csoilL_leaf=cpool_soil(1)
       csoilL_root=cpool_soil(2)
       !write(*,*) 'b', csoilL_leaf, csoilL_root
       csoilF_leaf=cpool_soil(3)
       csoilF_root=cpool_soil(4)

       csoilHs=cpool_soil(5)
       csoilHu=cpool_soil(6)
       csoilL=cpool_soil(7) 
       csoilF=cpool_soil(8)

       nsoilL_leaf=npool_soil(1)
       nsoilL_root=npool_soil(2)

       nsoilF_leaf=npool_soil(3)
       nsoilF_root=npool_soil(4) 
   
       nsoilHs=npool_soil(5)
       nsoilHu=npool_soil(6) 
       nsoilL=npool_soil(7)
       nsoilF=npool_soil(8) 

       litter_inp_leaf=0.0
       litter_inp_root=0.0
       n_cont_leaf0=0.0
       n_cont_root0=0.0
       ash_cont_leaf0=0.0
       ash_cont_root0=0.0


       do p=1, nv2
!     ab_litter(p)-percent of aboveground litter in total litter
!     n_littera(p)-N content in aboveground litter, %
!     ash_littera(p)-ash nontent in aboveground litter, % 
!       write(*,*)'VEGDAT',ab_litter(p),n_littera(p),n_litterb(p)
      !write(*,*) p, litter_in(p)
      litter_inp_leaf=litter_inp_leaf+litter_in(p)*ab_litter(p)
       litter_inp_root=litter_inp_root+litter_in(p)*(1-ab_litter(p))

      n_cont_leaf0=n_cont_leaf0+litter_in(p)*ab_litter(p)*n_littera(p)
      n_cont_root0=n_cont_root0+litter_in(p)*(1-ab_litter(p))*
     & n_litterb(p)
      ash_cont_leaf0=ash_cont_leaf0+litter_in(p)*ab_litter(p)*
     & ash_littera(p)
      ash_cont_root0=ash_cont_root0+litter_in(p)*(1-ab_litter(p))*
     & ash_litterb(p)
       end do

      if (litter_inp_leaf.ne.0.0)then
        n_cont_leaf0=n_cont_leaf0/litter_inp_leaf
      else
        n_cont_leaf0=0.0 
      end if
      if(litter_inp_root.ne.0.0)then
        n_cont_root0=n_cont_root0/litter_inp_root 
      else
        n_cont_root0=0.0
      end if
      if(litter_inp_leaf.ne.0.0)then
        ash_cont_leaf0=ash_cont_leaf0/litter_inp_leaf
      else
        ash_cont_leaf0=0.0
      end if
      if (litter_inp_root.ne.0.0)then
        ash_cont_root0=ash_cont_root0/litter_inp_root
      else
        ash_cont_root0=0.0
      end if
      if((litter_inp_leaf+litter_inp_root).ne.0.0)then
      n_cont_F0=(litter_inp_leaf*n_cont_leaf0+
     & litter_inp_root*n_cont_root0)/
     & (litter_inp_leaf+litter_inp_root)
      ash_cont_F0=(litter_inp_leaf*ash_cont_leaf0+
     & litter_inp_root*ash_cont_root0)/
     & (litter_inp_leaf+litter_inp_root)
      else
        n_cont_F0=0.0
        ash_cont_F0=0.0
      end if

       csoilL_leaf=csoilL_leaf+litter_inp_leaf*dt
      ! write(*,*)'m', litter_inp_leaf, csoilL_leaf, dt
      csoilL_root=csoilL_root+litter_inp_root*dt
      

       nsoilL_leaf=nsoilL_leaf+litter_inp_leaf*
     & n_cont_leaf0/(0.5*100)*dt
        nsoilL_root=nsoilL_root+litter_inp_root*
     & n_cont_root0/(0.5*100)*dt

! *0.5 to get per dry mass from per C
      if (csoilL_leaf.ne.0.0)then
        n_cont_leaf=100*0.5*nsoilL_leaf/csoilL_leaf
      else
        n_cont_leaf=0.0
      end if
      if (csoilL_root.ne.0.0)then
        n_cont_root=100*0.5*nsoilL_root/csoilL_root
      else
        n_cont_root=0.0
      end if
      if (csoilF.ne.0.0)then
        n_cont_F=100*0.5*nsoilF/csoilF
      else
        n_cont_F=0.0
      end if

      if (csoilF_leaf.ne.0.0)then
        n_cont_F_leaf=100*0.5*nsoilF_leaf/csoilF_leaf
      else
        n_cont_F_leaf=0.0
      end if
      if (csoilF_root.ne.0.0)then
        n_cont_F_root=100*0.5*nsoilF_root/csoilF_root
      else
        n_cont_F_root=0.0
      end if


       call koeffsom1_3a(n_cont_leaf0, ash_cont_leaf0, k1_som_leaf,
     & k3_som_leaf)
       call koeffsom1_3b(n_cont_root0,ash_cont_root0,k1_som_root,
     & k3_som_root)
       call koeffsom2a(n_cont_leaf0,ash_cont_leaf0,k2_som_leaf)
       call koeffsom2b(n_cont_root0,ash_cont_root0,k2_som_root)
       call koeffsom4_7(n_cont_F0,ash_cont_F0,k4_som,k5_som,
     &  k7_som) 


!******************************NO WOOD IN THIS VERSION
!c          if (temp.GT.0) then
!c                 k1_som_wood=0.0003
!c k3_som_wood=0.00008961
!c          k3_som_wood=0
!c                 
!c                        else
!c          k1_som_wood=0
!c                k3_som_wood=0
!c                k2_som_wood=0.000496
!c        end if
!* By now the temperature and moisture effect on wood decomposition are not considered
!* because fungi can regulate T and W themself

        call temdec1_3 (temp_lit,masswc_lit,k1_som_leaf,k3_som_leaf)  
        call temdec1_3 (temp_org,masswc_org,k1_som_root,k3_som_root)   
        call temdec2 (temp_org,masswc_org,k2_som_leaf)
        call temdec2 (temp_org,masswc_org,k2_som_root)
        call temdec4_7 (temp_org,masswc_org,k4_som, k7_som)
        call temdec6 (temp_min,masswc_min, k5_som, k6_som)
        call distrib_ab(n_cont_F0,k4_som,k4_som_a,k4_som_b)

!*Defines N saving coefficients
        Ml_leaf=0.1
        Ml_root=0.1

        call M_coeff(n_cont_leaf0, n_cont_F_leaf, Mf_leaf)
        call M_coeff(n_cont_root0, n_cont_F_root, Mf_root)

       
!c          call M_coeff(n_cont_F0, n_cont_F, Mf)

        if ((csoilHs/nsoilHs).LE.8.0) then
          Mhs=1.0
        else
          Mhs=0.8
        end if
        if ((csoilHu/nsoilHu).LE.8.0) then
          Mhu=1.0
        else
          Mhu=0.8
        end if


       csoilHs=csoilHs+((k4_som*bbact+k5_som*blumb)*nsoilF_root+
     & (k4_som_b*bbact+k5_som*blumb)*nsoilF_leaf-
     & k6_som*csoilHs)*deltat
     
!*/stable humus Mineral layer
       nsoilHs=nsoilHs+((k4_som+k5_som)*nsoilF_root*mt+
     & (k4_som_b+k5_som)*nsoilF_leaf*mt-k6_som*Mhs*nsoilHs)*deltat
 
       csoilHu=csoilHu+(k4_som_a*bbact*nsoilF_leaf-k7_som*csoilHu)*
     &  deltat
!*/unstable humus Organic layer
       nsoilHu=nsoilHu+(k4_som_a*nsoilF_leaf*mt-k7_som*Mhu*nsoilHu)*
     &  deltat
 
       csoilF_leaf=csoilF_leaf+((k3_som_leaf*csoilL_leaf)-
     &  (k2_som_leaf+k4_som+k5_som)*csoilF_leaf)*deltat

       nsoilF_leaf=nsoilF_leaf+((k3_som_leaf*nsoilL_leaf)-
     &    (k2_som_leaf*Mf_leaf+k4_som+k5_som)*nsoilF_leaf)*deltat

       csoilF_root=csoilF_root +((k3_som_root*csoilL_root)-
     &    (k2_som_root+k4_som+k5_som)*csoilF_root)*deltat

       nsoilF_root=nsoilF_root+((k3_som_root*nsoilL_root)-
     &    (k2_som_root*Mf_root+k4_som+k5_som)*nsoilF_root)*deltat
 
   

        csoilF=csoilF_leaf+csoilF_root
         nsoilF=nsoilF_leaf+nsoilF_root


        csoilL_leaf=csoilL_leaf-((k1_som_leaf+
     &   k3_som_leaf)*csoilL_leaf)*deltat

 
       
        csoilL_root=csoilL_root-((k1_som_root+
     &    k3_som_root)*csoilL_root)*deltat
 
        nsoilL_leaf=nsoilL_leaf-((k1_som_leaf*Ml_leaf+
     &    k3_som_leaf)*nsoilL_leaf)*deltat
      
  

        nsoilL_root=nsoilL_root-((k1_som_root*Ml_root+
     &    k3_som_root)*nsoilL_root)*deltat
             
        csoilL=csoilL_leaf+csoilL_root
        nsoilL=nsoilL_leaf+nsoilL_root
      
       if(csoilF.ne.0.0)then
       cflux=(k1_som_leaf*csoilL_leaf+k1_som_root*csoilL_root+
     &k2_som_leaf*csoilF_leaf+k2_som_root*csoilF_root+
     &k4_som*(csoilF_leaf-bbact*nsoilF_leaf)+
     &k4_som*(csoilF_root-bbact*nsoilF_root)+
     &k5_som*(csoilF_leaf-blumb*nsoilF_leaf)+
     &k5_som*(csoilF_root-blumb*nsoilF_root)+
     &k6_som*csoilHs+k7_som*csoilHu)*deltat/dt    !in kgCm-2s-1

        else
      cflux=(k1_som_leaf*csoilL_leaf+k1_som_root*csoilL_root+
     &k2_som_leaf*csoilF_leaf+k2_som_root*csoilF_root)*deltat/dt 
      end if
 
        !  write(*,*)'a',csoilL_leaf, csoilL_root
       cpool_soil(1)=csoilL_leaf
       cpool_soil(2)=csoilL_root


       cpool_soil(3)=csoilF_leaf
       cpool_soil(4)=csoilF_root

       cpool_soil(5)=csoilHs
       cpool_soil(6)=csoilHu
       cpool_soil(7)=csoilL 
       cpool_soil(8)=csoilF

       npool_soil(1)=nsoilL_leaf
       npool_soil(2)=nsoilL_root
    
       npool_soil(3)=nsoilF_leaf
       npool_soil(4)=nsoilF_root
      
         npool_soil(5)=nsoilHs
         npool_soil(6)=nsoilHu
       npool_soil(7)=nsoilL
       npool_soil(8)=nsoilF
        
          RETURN
        end  

!* Functions describing nitrogen saving as a function of initial N content in litter

           subroutine M_coeff (n_cont0, n_cont, M_c)
                 real n_cont, n_cont0, M_c
                 if (n_cont.LE.(1.16*n_cont0+0.44)) then
             M_c=0.1
                 else if (n_cont.GT.(1.16*n_cont0+0.44).AND.n_cont.LE.
     &     (1.16*n_cont0+1.5)) then
                 M_c=0.5
                 else
              M_c=1.0
             end if
           end
!*/////////////////////////////////// 

           subroutine  koeffsom1_3a( n_cont, ash_cont, k1_som, k3_som) 
           real n_cont, ash_cont, k1_som, k3_som
             k1_som=0.0005+0.0054*n_cont
             k3_som=0.0089+0.0078*n_cont
!c             k1_som=0.00466*n_cont
!c           k3_som=0.0041+0.0086*n_cont-0.00172*n_cont*n_cont
             if(k3_som.LT.0) then
                 k3_som=0
             end if
                  end



           subroutine  koeffsom1_3b(n_cont, ash_cont, k1_som, k3_som) 
             real n_cont, ash_cont, k1_som, k3_som
!c             k1_som=0.0005+0.0054*n_cont
!c             k3_som=0.0089+0.0078*n_cont
        !as aboveground
           k1_som=0.0136+0.0006*ash_cont
                 k3_som=0.0394-0.0021*ash_cont 
!c           k1_som=0.001795*n_cont+0.00177*n_cont*n_cont-0.0002222
!c             k3_som=0.0060168+0.000186*ash_cont+0.0057772*n_cont
           if(k1_som.LT.0) then
                 k1_som=0
                  end if
           end

           subroutine  koeffsom2a( n_cont, ash_cont, k2_som)
             real n_cont, ash_cont, k2_som 
!c             k2_som=0.000496
           k2_som=0.00060 
           end

           subroutine  koeffsom2b( n_cont, ash_cont, k2_som) 
           real n_cont, ash_cont, k2_som
!c             k2_som=0.0003578-0.000001617*ash_cont+0.00004894*n_cont
           k2_som=0.00126
           end

           subroutine  koeffsom4_7(n_cont, ash_cont, k4_som, k5_som, 
     &         k7_som) 
           real n_cont, ash_cont, k4_som, k5_som,k50_som, k7_som
               if(n_cont.LE.2.0) then
           k4_som=0.0005*n_cont
             else
           k4_som=0.001
             end if
             if(n_cont.LE.0.5)then
           k50_som=0.0
             else if (n_cont.GT.0.5.AND.n_cont.LE.2.0)then
           k50_som=2*0.007*n_cont/3-0.007/3
             else
           k50_som=0.007
             end if
             if(ash_cont.LE.5.0) then
           k5_som=0.2*ash_cont*k50_som
             else
           k5_som=k50_som
             end if
!*             if (100/n_cont.GT.35) then
!*             k5_som=0 !NO eaarthworms
!*             end if
             k7_som=0.00016
           end

           subroutine  distrib_ab(n_cont, k4_som, k4_som_a, k4_som_b)
             real n_cont, k4_som, k4_som_a, k4_som_b
             if(100*0.5/n_cont.LE.10) then
           k4_som_a=0
           k4_som_b=k4_som
            else if (100*0.5/n_cont.GT.10.AND.
     &       100*0.5/n_cont.LE.35) then
           k4_som_b=k4_som*(1.4-0.04*100*0.5/n_cont)
           k4_som_a=k4_som-k4_som_b
             else
           k4_som_b=0
           k4_som_a=k4_som
             end if
           end
!*////// Temperature f(T) and moisture g(W) effect on decomposition

           subroutine  temdec1_3 (T,W, k1_som, k3_som) 
           real T, W, k1_som, k3_som
           real f1,f3,g1,g3
           if(T.GT.(-5).AND.T.LE.1) then
             f1=0.1595+0.0319*T
                    else if (T.GT.1.AND.T.LE.35) then
             f1=0.1754*exp(0.0871*T)
             else if (T.GT.35.AND.T.LE.60) then
             f1=8.791-0.1465*T
             else
             f1=0.0
             end if
             if (T.GT.-3.AND.T.LE.7) then
!* the lowest is -3 in original ROMUL formulation
             f3=1.3
           else if (T.GT.7.AND.T.LE.60) then
           f3=60*1.3/53-1.3*T/53
           else
           f3=0.0
             end if
!c             if (W.GT.7.AND.W.LE.30) then
!c             g1=W/23-7/23
!c           else if (W.GT.30.AND.W.LE.300)then
!c           g1=1.0
!c           else if (W.GT.300.AND.W.LE.600)then
!* g1=2-2*W/600
!c             g1=1.0
!c           else
!c           g1=0.0
!c             end if
           if(W.GE.7.AND.W.LT.125)then
             g1=0.00000453*W**(2.5492)
             else if (W.GE.125.AND.W.LT.400)then
             g1=1.0
             else if (W.GE.400)then
!c             g1=1.0027*0.99157**(W-400) !as in the book
           g1=1.0027*0.9995**(W-400)
             else
             g1=0.0
             end if

             g3=g1
             k1_som=k1_som*f1*g1
           k3_som=k3_som*f3*g3
             end

           subroutine  temdec2 (T,W,k2_som)
           real T,W,k2_som
             real f2,g2
                 if (T.LE.-5.OR.T.GT.60)then
             f2=0.0
             else if(T.GT.(-5).AND.T.LE.1)then
             f2=0.1595+0.0319*T
           else if (T.GT.1.AND.T.LE.25)then
           f2=0.1754*exp(0.0871*T)
           else if (T.GT.25.AND.T.LE.35)then
           f2=1.534
           else if (T.GT.35.AND.T.LE.60)then
           f2=3.690-0.0615*T
             end if
!c                    if (W.GT.7.AND.W.LE.50)then
!c             g2=0.0233*W-0.163
!c           else if (W.GT.50.AND.W.LE.90) then
!c           g2=1.312-0.00625*W
!c           else if (W.GT.90.AND.W.LE.1200)then
!c           g2=0.811-0.00068*W
!c           else
!c           g2=0.0
!c             end if

             if(W.GE.7.AND.W.LT.125)then
             g2=0.00000453*W**(2.5492)
             else if (W.GE.125.AND.W.LT.400)then
             g2=1.0
             else if (W.GE.400)then
!c             g1=1.0027*0.99157**(W-400) !as in the book
           g2=1.0027*0.998**(W-400)
             else
             g2=0.0
             end if


             k2_som=k2_som*f2*g2
           end

                         
          subroutine  temdec4_7 (T, W, k4_som,k7_som) 
             real T, W, k4_som,k7_som
             real f7,f4,g7,g4
          
             if (T.LE.-5.OR.T.GT.60)then
             f7=0.0
             else if(T.GT.(-5).AND.T.LE.1)then
             f7=0.1595+0.0319*T
           else if (T.GT.1.AND.T.LE.25)then
           f7=0.1754*exp(0.0871*T)
           else if (T.GT.25.AND.T.LE.35)then
           f7=1.534
           else if (T.GT.35.AND.T.LE.60)then
           f7=3.690-0.0615*T
             end if
             if(T.GT.-5.AND.T.LE.1)then
                 f4=0.1595+0.0319*T 
           else if (T.GT.1.AND.T.LE.20) then
           f4=0.1754*exp(0.0871*T)
           else if (T.GT.20.AND.T.LE.40)then
           f4=1
             else if (T.GT.40.AND.T.LE.80)then
           f4=2.0-0.025*T
             else 
           f4=0.0
             end if

        
             if (W.GT.7.AND.W.LE.50)then
             g7=0.0233*W-0.163
           else if (W.GT.50.AND.W.LE.90)then
           g7=1.312-0.00625*W
           else if (W.GT.90.AND.W.LE.1200)then
           g7=0.811-0.00068*W
           else
           g7=0.0
             end if
             if (W.LE.40)then
             g4=0.025*W
           else if (W.GT.40.AND.W.LE.400)then
           g4=1.0
           else if (W.GT.400.AND.W.LE.700)then
           g4=2.333-0.0033*W
           else
           g4=0.0
             end if
          
                 k4_som=k4_som*f4*g4
             k7_som=k7_som*f7*g7
             end

           subroutine  temdec6 ( T,W, k5_som,k6_som) 
           real T,W, k6_som, k5_som
             real f6,g6,f5, g5


             if (T.LE.-5.OR.T.GT.50)then
             f5=0.0
             else if(T.GT.(-5).AND.T.LE.1)then
             f5=0.078+0.0156*T 
           else if (T.GT.1.AND.T.LE.13) then
             f5=0.0675*exp(0.2088*T)
           else if (T.GT.13.AND.T.LE.25) then
           f5=1.0
           else if (T.GT.25.AND.T.LE.50)then
           f5=2.0-0.04*T
           end if
             if (W.GT.2.AND.W.LE.15) then
             g5=0.0769*W-0.1538
           else if (W.GT.15.and.W.LE.70) then
           g5=1.0
           else if (W.GT.70.AND.W.LE.120) then
           g5=2.4-0.02*W
           else
           g5=0.0
             end if
             if (T.LE.-5.OR.T.GT.60) then
             f6=0.0
             else if(T.GT.(-5).AND.T.LE.1)then
             f6=0.1595+0.0319*T
           else if (T.GT.1.AND.T.LE.27.5)then
           f6=0.1754*exp(0.0871*T)
           else if (T.GT.27.5.AND.T.LE.35)then
           f6=1.95
           else if (T.GT.35.AND.T.LE.60)then
           f6=4.68-0.078*T
           end if
             if (W.LE.40) then
             g6=0.025*W
             else
           g6=1 
             end if
             k5_som=k5_som*f5*g5
                  k6_som=k6_som*f6*g6
             end
        
